Hands-on Exercise 7c: Analytical Mapping

Author

Goh Si Hui

Published

February 26, 2024

Modified

February 26, 2024

1 About this Exercise

In this exercise, we will learn how to plot analytical maps.

2 Getting Started

2.1 Installing and Loading R Packages

For this exercise, other than tmap, we will use the following packages:

  • tidyverse for tidying and wrangling data
  • sf for handling geospatial data

The code chunk below uses p_load() of pacman package to check if the abovementioned packages are installed in the computer. If they are, they will be launched in R. Otherwise, pacman will install the relevant packages before launching them.

Show the code
pacman::p_load(tidyverse, tmap, sf)

2.2 Importing the Data

For the purpose of this hands-on exercise, a prepared data set called NGA_wp.rds will be used. The data set is a polygon feature data.frame providing information on water point of Nigeria at the LGA level.

We will use read_rds() function to import the data into R.

Show the code
nga <- read_rds("data/rds/NGA_wp.rds")
glimpse(nga)
Rows: 774
Columns: 9
$ ADM2_EN          <chr> "Aba North", "Aba South", "Abadam", "Abaji", "Abak", …
$ ADM2_PCODE       <chr> "NG001001", "NG001002", "NG008001", "NG015001", "NG00…
$ ADM1_EN          <chr> "Abia", "Abia", "Borno", "Federal Capital Territory",…
$ ADM1_PCODE       <chr> "NG001", "NG001", "NG008", "NG015", "NG003", "NG011",…
$ geometry         <MULTIPOLYGON [m]> MULTIPOLYGON (((548795.5 11..., MULTIPOL…
$ total_wp         <int> 17, 71, 0, 57, 48, 233, 34, 119, 152, 66, 39, 135, 63…
$ wp_functional    <int> 7, 29, 0, 23, 23, 82, 16, 72, 79, 18, 25, 54, 28, 55,…
$ wp_nonfunctional <int> 9, 35, 0, 34, 25, 42, 15, 33, 62, 26, 13, 73, 35, 36,…
$ wp_unknown       <int> 1, 7, 0, 0, 0, 109, 3, 14, 11, 22, 1, 8, 0, 37, 88, 1…

3 Basic Choropleth Maps

We will plot

total <- tm_shape(nga) +
  tm_polygons("total_wp",
              palette = "Greens",
              style = "equal",
              lwd = 0.1,
              alpha = 1) +
  tm_layout(main.title = "Distribution of All Water Points",
            main.title.size = 1,
            legend.outside = FALSE)

total

functional <- tm_shape(nga) +
  tm_polygons("wp_functional",
              palette = "Blues",
              style = "equal",
              lwd = 0.1,
              alpha = 1) +
  tm_layout(main.title = "Distribution of Functional Water Points",
            main.title.size = 1,
            legend.outside = FALSE)

functional

tmap_arrange(total, functional, nrow = 1)

4 Choropleth Map for Rates

4.1 Deriving Proportion of Functional Water Points and Non-Functional Water Points

nga <- nga %>%
  mutate(pct_functional = round(wp_functional/total_wp,2)) %>%
  mutate(pct_nonfunctional = round(wp_nonfunctional/total_wp,2))

glimpse(nga)
Rows: 774
Columns: 11
$ ADM2_EN           <chr> "Aba North", "Aba South", "Abadam", "Abaji", "Abak",…
$ ADM2_PCODE        <chr> "NG001001", "NG001002", "NG008001", "NG015001", "NG0…
$ ADM1_EN           <chr> "Abia", "Abia", "Borno", "Federal Capital Territory"…
$ ADM1_PCODE        <chr> "NG001", "NG001", "NG008", "NG015", "NG003", "NG011"…
$ geometry          <MULTIPOLYGON [m]> MULTIPOLYGON (((548795.5 11..., MULTIPO…
$ total_wp          <int> 17, 71, 0, 57, 48, 233, 34, 119, 152, 66, 39, 135, 6…
$ wp_functional     <int> 7, 29, 0, 23, 23, 82, 16, 72, 79, 18, 25, 54, 28, 55…
$ wp_nonfunctional  <int> 9, 35, 0, 34, 25, 42, 15, 33, 62, 26, 13, 73, 35, 36…
$ wp_unknown        <int> 1, 7, 0, 0, 0, 109, 3, 14, 11, 22, 1, 8, 0, 37, 88, …
$ pct_functional    <dbl> 0.41, 0.41, NaN, 0.40, 0.48, 0.35, 0.47, 0.61, 0.52,…
$ pct_nonfunctional <dbl> 0.53, 0.49, NaN, 0.60, 0.52, 0.18, 0.44, 0.28, 0.41,…
functional_rate <- tm_shape(nga) +
  tm_polygons("pct_functional",
              n = 10,
              palette = "Blues",
              style = "equal",
              lwd = 0.1,
              alpha = 1,
              legend.hist = TRUE) +
  tm_layout(main.title = "Rate Map of Functional Water Points by LGAs",
            main.title.size = 1,
            legend.outside = TRUE)

functional_rate

5 Extreme Value Maps

5.1 Percentile Map

5.1.1 Data Preparation

nga <- nga %>%
  drop_na()
percent <- c(0, 0.01, 0.1, 0.5, 0.9, 0.99, 1)
var <- nga["pct_functional"] %>%
  st_set_geometry(NULL)

quantile(var[,1], percent)
  0%   1%  10%  50%  90%  99% 100% 
0.00 0.00 0.22 0.48 0.86 1.00 1.00 

5.2 Creating get.var function

get.var <- function(vname, df) {
  v <- df[vname] %>%
    st_set_geometry(NULL)
  v <- unname(v[,1])
  return(v)
}

5.3 Percentile Mapping Function

percentmap <- function(vnam, df, legtitle=NA, mtitle="Percentile Map"){
  percent <- c(0,.01,.1,.5,.9,.99,1)
  var <- get.var(vnam, df)
  bperc <- quantile(var, percent)
  tm_shape(df) +
  tm_polygons() +
  tm_shape(df) +
     tm_fill(vnam,
             title=legtitle,
             breaks=bperc,
             palette="Blues",
          labels=c("< 1%", "1% - 10%", "10% - 50%", "50% - 90%", "90% - 99%", "> 99%"))  +
  tm_borders() +
  tm_layout(main.title = mtitle, 
            title.position = c("right","bottom"))
}
percentmap("total_wp", nga)

6 Box Map

ggplot(data = nga,
       aes(x = "",
           y = wp_nonfunctional)) +
  geom_boxplot()

boxbreaks <- function(v,mult=1.5) {
  qv <- unname(quantile(v))
  iqr <- qv[4] - qv[2]
  upfence <- qv[4] + mult * iqr
  lofence <- qv[2] - mult * iqr
  # initialize break points vector
  bb <- vector(mode="numeric",length=7)
  # logic for lower and upper fences
  if (lofence < qv[1]) {  # no lower outliers
    bb[1] <- lofence
    bb[2] <- floor(qv[1])
  } else {
    bb[2] <- lofence
    bb[1] <- qv[1]
  }
  if (upfence > qv[5]) { # no upper outliers
    bb[7] <- upfence
    bb[6] <- ceiling(qv[5])
  } else {
    bb[6] <- upfence
    bb[7] <- qv[5]
  }
  bb[3:5] <- qv[2:4]
  return(bb)
}
get.var <- function(vname,df) {
  v <- df[vname] %>% st_set_geometry(NULL)
  v <- unname(v[,1])
  return(v)
}
var <- get.var("wp_nonfunctional", nga) 
boxbreaks(var)
[1] -56.5   0.0  14.0  34.0  61.0 131.5 278.0
boxmap <- function(vnam, df, 
                   legtitle=NA,
                   mtitle="Box Map",
                   mult=1.5){
  var <- get.var(vnam,df)
  bb <- boxbreaks(var)
  tm_shape(df) +
    tm_polygons() +
  tm_shape(df) +
     tm_fill(vnam,title=legtitle,
             breaks=bb,
             palette="Blues",
          labels = c("lower outlier", 
                     "< 25%", 
                     "25% - 50%", 
                     "50% - 75%",
                     "> 75%", 
                     "upper outlier"))  +
  tm_borders() +
  tm_layout(main.title = mtitle, 
            title.position = c("left",
                               "top"))
}
tmap_mode("view")
boxmap("wp_nonfunctional", nga)
tmap_mode("plot")